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Fourier-Legendre decomposition (FLD) of solar Doppler imaging data is a promising method to estimate the sub-surface 
solar meridional flow. FLD is sensible to low-degree oscillation modes and thus has the potential to probe the deep 
meridional flow. We present a newly developed code to be used for large scale FLD analysis of helioseismic data as 
provided by the Global Oscillation Network Group (GONG), the Michelson Doppler Imager (MDI) instrument, and the 
upcoming Helioseismic and Magnetic Imager (HMI) instrument. First results obtained with the new code are qualitatively 
comparable to those obtained from ring-diagram analyis of the same time series. 
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1 Introduction 

From summer 2010 on, the Helioseismic and Magnetic Im- 
ager (HMI) instrument aboard the Solar Dynamics Observa- 
tory (SDO) satellite is expected to deliver continuous high 
resolution (4 k x 4 k) full-disk Doppler images at a cadence 
of 45 seconds. The availability of this data will surely mark 
a new era for helioseismology, but the huge amount of data 
that HMI will produce also poses a challenge for helioseis- 
mic data processing pipelines. With the given resolution and 
cadence, the HMI full-disk Doppler images alone will sum 
up to roughly 2 terabytes per month. 

The current major sources of solar full-disk Doppler 
data are the Globa l Oscillation Network Group (GONG; 
Harvey etal.1 1996|) and t he Michelson Doppler Imager 



ey < 



(MDI; Scherrer et al. 1995) instrument aboard the Solar and 



Heliospheric Observatory (SOHO) spacecraft. They deliver 
Dopplergrams with a maximum volume of about 100 giga- 
bytes per month, each. Both GONG and MDI have been 
operating for more than one decade and it is exciting that 
one can now apply helioseismic techniques to data spanning 
about one solar cycle from two independent observatories. 

The characteristics of the solar meridional circulation 
are an important ingre dient to some solar dynam o models, 
see e.g. the review by Dikpati & Gilmanl ( 20091) and refer- 
ences therein. A poleward flow of the order of 10 to 20 m/s 
can be derived at the surface and in the near-surface layers 
by a variety of techniques (e.g. Duvalll 1979 : Haber et al 



20021: lHathawavlll996t iKomm et alJll993l IWohl & Braisa 



2001) but no significant evidence for a return flow in deeper 
layers has been detected so far. 

Fourier-Legendre spectral decomposition (FLD) is sen- 
sitive to low-degree modes and thus has the potential to de- 



tect flows in d eep layers. From encouraging results of pre- 
vious research ( Braun & Fan! 1998 : Krieger et al.ll2007 ) and 
the progress of our ongoing work, we are confident that the 
Fourier-Legendre analysis can become a standard tool for 
meridional flow measu rements be low the surface along with 
ring-diagra m analysis (Hill 19881) and time-distance helio- 
seismology dDuvall et al l ll993l) . 

However, this requires a data processing pipeline that is 
able to handle efficiently the available and upcoming helio- 
seismic data. In this paper, we present a newly developed 
code for Fourier-Legendre analysis that fits these require- 
ments. The code has been subject to several tests. Here we 
present first results obtained with this code and compare 
them to results from ring-diagram analysis of the same time 
series. It is our intention making the code available to the 
helioseismology community as part of the European Helio- 
and Asteroseismology Network (HELAS). 

2 Methods 

2.1 The Fourier-Legendre spectral decomposition 

The "Fourier-Hankel spectral decomposition" technique 
has been used to s t udy jo-mode scatt e ring b y sunspots 



(IBogdan et al.1 Il993t iBraun et al.l 119871 1 19921) . Later the 



method was also applied for sub-surfa ce m eridional flow 
measu rements by Braun & Fan! ( 1998 ) and Krieger et al. 
d2007l) . 



In plane geometry, Hankel functions are a good approx- 
imation to Legendre functions and in some of the publica- 
tions mentioned above they were used instead of Legendre 
functions because they are numerically much easier to com- 
pute. In this work we do not use Hankel functions, and to 
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avoid confusion we prefer to use the more descriptive no- 
tation "Fourier-Legendre decomposition" for the remainder 

of this paper. 

Following iBraunetaD dl988l) and lBraun & Fan! dl998l) . 
the surface oscillation signal SV(6, <fi, t) within an annular 
region around a point of interest can be represented as a 
superposition of inward and outward traveling waves of the 
form 



SV(6, & t) = J2 [Ai mv xr(6) 



Imv 



+ B lmu (XJ n )* (9)] e^™^ 2 ™*), (1) 

where Ai mv and Bi mv are the complex amplitudes of the 
oscillation modes of the temporal frequency v, the harmonic 
degree I and azimuthal order m. The angle 9 equals zero in 
the center of the region of interest and increases with radial 
distance from the center. The basis functions 



XP(8) = Nf 



Pr{cos9) - —Qr(cost 

7T 



(2) 



and its complex conjugate {Xf 1 )* are superpositions of the 
associated Legendre functions of first and second kind P" 1 



and Q?\ and Nf 



( i\m (l-m)\ 
\ l ) (l+mS\ 



is a normalization factor. 



For meridional flow measurements, the center of the an- 
nular region is identical to either the northern or southern 
pole, so that Ai„ w and Bi mv become the poleward and equ- 
atorward waves, respectively. 

Using equation (Q~|), the mode amplitudes can be extract- 
ed from the measured Doppler signal by 

T 2tt max 
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x X z m (0)0d0#di, (4) 

where C; is a normalization factor which is approximately 
given byny/l(l + l)/[2(0 max -0 min )], max and min mark 
the latitudinal extent of the annular region, and T is the total 
length of the observed time series. 

3 The numerical code 

The decomposition code is written in the C language for ef- 
ficiency and portability. It is parallelized using the Message 
Passing Interface (MPI) standard. Further steps in the post- 
processing pipeline such as peak fitting and the inversion 
are currently implemented as prototype code in IDL (Inter- 
active Data Language). 



The design goals for the code were: 1) to address the 
most relevant issues and open questions that arise from pre- 
vious publications, 2) to provide the flexibility to easily use 
data from common sources such as MDI, GONG and the 
upcoming HMI, 3) to write a code that is fast enough to 
process all available data in reasonable time, and 4) to make 
the code flexible enough to not only use it for meridional 
flow measurements, but also to measure jo-mode ab sorption 
in sunspots as demonstrated bv lBraun et al.l (119881). 



We ad dressed the issues as discussed by IBraun & Fan 



(119981) and lKrieger et al.l (120071) . This includes proper han- 
dling of the spherical geometry by using Legendre func- 
tions instead of Hankel functions (this was already done by 
Braun & Fanl though) and an inversion for the flow-profile. 
Requirement 2) is implemented by following a modular 
concept: for each data source, a plug-in module is provided 
that knows how to read data from the respective source and 
forwards the Doppler images to the decomposition module 
in a standardized format. Goal 4) will be implemented in the 
future. 

On a modern quad-core processor, one month of GONG 
data can be processed in less than one day with a setup as we 
use it in this paper (see below). The decomposition into the 
mode coefficients only takes a minor fraction of the total 
runtime and the further steps are not optimized for speed 
yet, so that we expect future versions of the code to run 
much faster. The code also scales well up to the amounts of 
data we expect from the HMI instrument later this year. 

3.1 Data preparation 

The Doppler images are processed in chunks. Each chunk 
consists of a configurable number of Dopplergrams which 
usually equals to 60 for the one-minute cadence of GONG 
and MDI. 

The Dopplergrams are interpolated to an equidistant he- 
liocentric 6, 4> -grid. Effects of the solar rotation and rela- 
tive movement between instrument and Sun are removed by 
subtracting the mean Dopplergram of one chunk from each 
Dopplergram in the chunk. Each Dopplergram is then apod- 
ized using a Hann window function to avoid spatial alias- 
ing in the resulting power spectra. Bad or missing Doppler- 
grams are detected and the mode amplitudes of the corre- 
sponding timestep are set to zero. A binary mask is stored 
along with the time series that allows to properly detect and 
account for gaps in the time series. 

3.2 Decomposition & frequency fitting 

For the decomposition of the mode co efficien ts, we fol- 



low an approach that was first used by iBrownl d 19851) for 
fast Spherical Harmonic decomposition of solar oscillation 
data. The integration in longitude in Eqs. (0 and (HJi is im- 
plemented as a Fast Fourier Transform (FFT). Compared 
to a straightforward numerical evaluation of the integrals, 
this results in a greatly improved performance of the overall 
computation. 
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Because the associated Legendre functions P[ n and Q™ 
can only be calculated reliably using a recursion relation 



Table 1 Latitudinal position and width of the patches in 
heliocentric coordinates 



( ZhangHl996h . we precompute the basis functions X™ for 



all combinations of I, m and 0. Since we are only interested 
in low azimuthal orders m = —25 . . .+25, the lookup tables 
consume a few hundred megabytes and fit easily into main 
memory. 

The resulting time series for the complex amplitudes 
Ai m and Bi m are stored as FITS (Flexible Image Trans- 
port System) binary tables for the further steps in the pro- 
cessing pipeline. From the power spectra of Ai mv and Bi mv 
the peak frequencies are determined by fitting asymmetric 
Lorentzian profiles. From these fits the frequency differ- 
ences between the poleward and equatorward propagating 
waves are determined. 

3.3 Inversion 

The frequency shift between the poleward- and equator- 
ward propagating waves is used to invert for the hori- 
zontal component of the sub-surface meridional flow. The 
frequency shift Ay is r elated to the meridional flow by 
dGough&Toomrel 19831) 



I J(U mcT (r))K nl {r) dr 
_o 

oo 

ttRq / K n i (r) dr 
o 



(5) 



where R@ is the solar radius, (U meY (r)) is the mean merid- 
ional flow averaged over the patch, and K n i(r) are the 
energy densi ty kernels, which are calculated from solar 
'Model S' bv lChristensen-Dalsgaard et al. dl996l) . Based on 
this equation, inversions for the meridional flow are carried 
out by employing a SOLA technique ( Subtractive Optimally 
Localized A veraging) as described in iPijpers & Thompson 
dl992l[l994l) . Taking measurement errors of the frequencies 
into account, the inverted solution has to be regularized. For 
simplicity, we used one regularization parameter for all po- 
sitions in the Sun. We selected a rather large regularization 
parameter in order to give trust only to those modes where 
a precise frequency measurement was possible. 



center latitude [deg] 


longitudinal extent [deg] 


±45.0 


-30.0 to +30.0 


±37.5 


-45.0 to +45.0 


± 30.0 


-45.0 to +45.0 


± 22.5 


-52.5 to +52.5 


± 15.0 


-52.5 to +52.5 


±7.5 


-52.5 to +52.5 





-52.5 to +52.5 



These are the same parameters that were used for position- 
ing 16x16 degree square patches for the ring-diagram anal- 
ysis. For each patch the Fourier-Legendre decomposition is 
carried out and the meridional flow is determined from the 
frequency shifts between pole- and equatorward propagat- 
ing waves of harmonic degrees I = 100 - 1000 and radial 
order n — - 11. 

4.1 Comparison with ring-diagram analysis 

Figure Q] shows cross sections of the flow profiles at depths 
of 3, 5 and 7 Mm obtained from Fourier-Legendre decom- 
position (left panel) and ring-diagram analysis (right panel). 
With both techniques the measurement error is lowest close 
to the surface and at the equator (±0.2 m/s for FLD and 
ring-diagram analysis) and highest in the deep layers and 
at high latitudes (±0.5 m/s for ring-diagram analysis and 
0.8 m/s for FLD). 

Both methods reproduce the same qualitative features 
of the flow profile. The direction of the flow is mainly pole- 
ward in the order of 20 m/s. FLD and ring-diagram both fa- 
vor a weak equator-crossing flow but with opposite sign. 
With both techniques the derived flow velocities increase 
with latitude. There is a clear asymmetry between both 
hemispheres visible in the velocity profiles from both meth- 
ods. On the northern hemisphere the velocities tend to de- 
crease with depth, whereas the curves for the different depth 
agree within their error marings down to -22.5 degree lati- 
tude. Interstingly, both methods show the same bend in the 
velocity profile at -7.5 degree latitude. 



4 Results 

For a first test, we use GONG data from January and Febru- 
ary 2006. These data sets where chosen because at the be- 
ginning of this period the GONG duty cycle was compa- 
rably high (91%) and it is in the middle of the declining 
phase of the past solar cycle, so that a minimum of tempo- 
ral change in the flow pattern can be estimated. 

To estimate the meridional flow as a function of depth 
and latitude, the Dopplergrams are subdivided into smaller 
patches. All patches have a height of 16 degree in latitude 
while the longitudinal extent depends on the position on 
the disk in order to avoid foreshortened regions. The loca- 
tion and width of the patches we used is shown in Table Q] 



5 Conclusions & Discussion 

Combined with helioseismic inversion techniques Fourier- 
Legendre decomposition of Doppler imaging data can be 
used to derive the sub-surface solar meridional flow. The 
method is sensible to low-degree modes and thus could 
greatly increase the range in depth that is currently acces- 
sible to other methods such as ring-diagram analysis. 

We developed a new numerical code suitable for 
Fourier-Legendre analysis of large sets of Dopplergrams as 
provided by GONG, MDI and the upcoming HMI instru- 
ment. For a first test of the new code, we compared the 
near-surface flow velocities derived from FLD to those ob- 
tained from ring-diagram analysis of the same data sets. 
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Fourier— Legendre Decomposition 



Ring Diagram Analysis 
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Figure 1 Near surface meridional flow as derived from Fourier-Legendre decomposition and ring-diagram analysis of 
two month of GONG data. Errorbars span the 3 a interval at each point. The graphs for 5 and 7 Mm are shifted by 0.5 and 
1.0 degrees respectively for better visibility. 



Both methods result in qualitatively comparable flow veloc- 
ities with small errors in the order of one percent. However, 
the absolute values for the velocities obtained with both 
methods are not in complete agreement. As for the inver- 
sion of the Fourier-Legendre data one regularization param- 
eter was used for all positions in the Sun, the discrepancies 
between the two methods might be due to this systematic 
effect. 

A careful analysis of possible sources of further sys- 
tematic errors needs to be carried out in future. For the 
case of FLD this also includes the impact of leakage be- 
tween modes of neighboring degree I and order m which 
is not taken into account in the present paper. Possible so- 
lutions are either to use only a limited set of values of I 
for whi ch orthogonality o f the Legendre functions is guar- 
anteed ( Braun et alJI 1988 ). Another possibility is to include 
the covariance matrix of the mode coefficients in the proce- 
dure for the frequency determination of the modes. 

Inversions for the meridional flow in deeper layers of 
the Sun as well as long-term studies of the variability of the 
flow will be carried out when the open issues are resolved 
in the near future. 
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